library(reshape2)
library(ggplot2)
theme_set(theme_classic() +
	theme(text=element_text('times')) +
	theme(strip.background=element_blank()) +
	theme(strip.text.x=element_text(face='bold',size=12)))

source('funk1.R')
# R1: 0   -(b)-> T
# R2: T   -(d)-> 0
# R3: T+V -(e)-> Y+V
# R4: Y   -(u)-> p*V
# R5: V   -(c)-> 0

params <- list(
	b=1, d=0.1, e=0.001, u=0.5, p=1000, c=0.5
	#b=1, d=0.1, e=0.001, u=0.5, p=1000, c=10
)
X0 <- c(T=params$b/params$d, Y=5, V=10)
#X0 <- c(T=params$b/params$d, Y=0, V=10)
nsim <- 20

X.ode <- funk1(X0, params, model.euler, 0.01)
X <- cbind(melt(as.data.frame(X.ode), id='time'), type='ode', group=0)

library(parallel)
cl <- makeCluster(detectCores())
clusterExport(cl=cl, varlist=c("funk1", "X0", "params", "model.solve", "model.tau_leaping", "model.maruyama", "model.euler"), envir=environment())
#X.ssa <- parLapply(cl, 1:nsim, function(i) funk1(X0, params, model.maruyama, 10))
X.ssa <- parLapply(cl, 1:nsim, function(i) funk1(X0, params, model.tau_leaping, 0.01))
stopCluster(cl)

for(i in 1:length(X.ssa))
	X <- rbind(X, cbind(melt(as.data.frame(X.ssa[[i]]), id='time'), type='ssa', group=i))

X <- rbind(X, cbind(aggregate(value ~ time+variable, subset(X,type=='ssa'), mean), type='avg', group=length(X.ssa)+1))

X$variable <- factor(X$variable,
	levels=c('T','Y','V'),
	labels=c('Target cells','Infected cells','Free virus'))

library(plyr)  # for subset, to limit yaxis

print(ggplot(X, aes(time,value,group=group,alpha=type,size=type,color=type)) +
	geom_line(subset=.(variable=='Target cells'|(value<5&variable=='Infected cells')|(value<15000&variable=='Free virus'))) +
	#geom_line() +
	facet_wrap(~ variable, scales='free_y') +
	expand_limits(y=0) +
	scale_alpha_manual(values=c(1,min(1,2.5/nsim),1)) +
	scale_size_manual(values=c(0.5,0.3,0.3)) +
	scale_color_manual(values=c('black','black','red')) +
	xlab('Time units') +
	ylab('Cell/viral units') +
	theme(legend.position='none'))
